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I. INTRODUCTION 

The prominent light-matter interaction of graphene 
has attracted immense interest due to the tunable plas¬ 
monic excitations in the THz^“® and mid-infrared^’® 
regimes. These charge density excitations in graphene, 
which are mostly explained by the density-density and 
current-current correlation functions of its tt electron 
gas®, are well-explored in both theoreticaland 
experimentah^^^® directions. 

The manipulation of the density of states (DOS) of the 
TT and TT* bands of graphene can be a tool for tailoring 
its plasmonic excitations—a scenario realizable through 
the exposure of graphene to either mechanical stress^^’^® 
or perpendicular magnetic field^®^®^. Also, as implied 
by the Pauli exclusion principle, the manipulation of the 
electronic occupation within the tt and tt* bands alters 
the response to the electromagnetic (EM) perturbations. 
Gate-controlled optical absorption of graphene^® and the 
broadband optical gain resulting from the inversion of 
the electronic population under femtosecond laser pulse 
irradiation®®^®® are examples of altering the EM response 
of graphene via pushing its electronic occupation into 
steady and transient nonequilibrium states, respectively. 

Modifying light absorption by electrical signals would 
integrate optics and electronics, a long-sought goal in 
plasmonics®^. Breaking the spatial-temporal symmetries 
would also open up the possibility of rectifying the plas¬ 
monic current to convert light signals into directed elec¬ 
trical signals.®®’®® The directional symmetry is most di¬ 
rectly broken by applying an electric field within the two- 
dimensional (2D) layer. This will modify the spectrum, 
alter the response of the system and induce non-linear 
and thermal effects, while the presence of electrical con¬ 
tacts can lead to Dyakonov-Shur instabilities®®^®®. Also, 
population inversion induced by optical pumping can 
lead to a negative total dynamic conductivity in graphene 
at THz/far-infrared frequencies paving the way towards 
graphene-based laser devices.®®’®^ 

In this work, we will investigate the interplay between 
the electrical conductivity and the plasmonic excitations 
in graphene samples by assuming a moderate electric flux 


across the sample. To present the essence of our work, we 
focus on the analysis of the linear, intravalley response of 
drifting tt electron gas to longitudinal EM perturbations, 

E;(r,t) = £;oe*(‘^-"--*) ■,Eo\\q (1) 

The paper is organized as follows. Sec. II contains the 
basics of the linear response theory of Dirac systems and 
its generalization to nonequilibrium systems. In Sec. Ill, 
we present the analytical approximation valid for small 
drift velocities at finite doping. In Sec. IV, a general 
discussion for doped systems is given, and in Sec. V, we 
present the numerical results for the case of zero doping. 
We close with a summary and conclusions. The paper is 
supplemented by four appendices which provide details 
on the analytical calculations. 


II. LINEAR RESPONSE OF A DRIVEN DIRAC 
SYSTEM 


Within the random phase approximation®® (RPA), the 
response of the tt electrons at equilibrium to longitudinal 
EM perturbations is mainly described by the intravalley 
dynamical polarization function of graphene^®’®®, 

n(g,w) = X / d^k{fs,s'{k,q) 

sp=± ^ ^ 2 ) 

{k + q)]-npiE’^jk)] 
E^'{k-\-q) — E’^ik) — huj — i0+ 


where gs{gv) = 2 denotes the spin (valley) degeneracy, the 
prefactor /s,s'(fc, q) represents the band overlap integral 
and E^{k) describes the energy dispersion of the valence 
(s = — 1) and conduction (s = 1) bands. Employing the 
tight-binding model, if accompanied by the Dirac cone 
approximation, yields A®(fc) = s{3ato/2)k together with. 


fs,s'{k, q) 
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2^ |fc + q| ^ 
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where a ~ 0.142nm, to ~ 2.7eV and k are respectively 
the carbon-carbon bond length, the nearest-neighbor 
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hopping amplitude and the crystal momentum measured 
with respect to the Dirac points. In addition, 0fc(q) is the 
angle between k (q) and e^- The equilibrium electronic 
occupation is described by the Fermi-Dirac statistics, i.e.. 


nplE] 


1 + exp 


fE — Ep 
^ kpE 


) 


-1 


(4) 


where Ep is the Fermi energy measured with respect to 
the neutrality point. The highest occupied eigen-states 
in the reciprocal space are located at circles centered at 
the Dirac points. The disk enclosed by such a circle is 
referred to as the Fermi disk with the Fermi wavevector, 
kp = 2 \Ep\ /{Sato), being its radius. 

In the presence of drift, the tt electron gas reaches a 
new equilibrium through gaining momentum and kinetic 
energy from the drain-source voltage, Vds , and losing part 
of it via electron scattering mechanisms^®’^^. Instead of 
finding the eigen-states and energy eigen-values of the 
Hamiltonian that includes (the local electric field 
corresponding to Vds) and the sources of scattering, we 
adopt a semi-classical approach in which the electronic 
occupation is altered by the drain-source voltage, while 
the crystal Hamiltonian and consequently its DOS, re¬ 
main intact. As a result of this approach, the drift- 
induced modification to the dynamical polarization of 
graphene, which is defined as the difference between the 
dynamical polarization in the presence of drift n*(q,w) 
by its no-drift counterpart n°(q,u;), is given by^®. 


J (fk{fs,s'{k,q) 

, (5) 

AnpjE^ (fc + g)] - Anp[E^{k)] 

E^'{k + q) — E^{k) — hjjj — *0+ 

where Anp[E] = n*p[E\—n°p[E\ denotes the drift-induced 
modification to the electronic occupation. 


III. ANALYTIC APPROXIMATION 


case of pure electron or hole transport and relegate the 
special case of doping-levels close to half-filling to Sec. 
V. Within the low-temperature {kpT <C \Ep\) and low- 
drift {kshift^kp) regime, which is a relevant regime for 
usual doping levels and current densities, the occupation 
function of the drifting electron gas rip [E] can be approx¬ 
imated via feeding Ep = Ep{kp/kp) from Eq. (6) into 
the Fermi-Dirac occupation function. This yields, 

Anp{E,k) ^-Ep[’^^]6{E - Ep) cosOk (7) 

kp 

with 5 being the Dirac delta function. Feeding Eq. (5) 
with this spike-like Anp{E,k), if accompanied by the 
Dirac cone approximation, yields an analytic expression 
for An(g,a;) (see Appendix A). For brevity, we present 
this analytic expression in terms of the dimensionless 
variables q = q/kp, uj = huj/Ep and uj' = uj + lO'*', 


An(g,w) 


D{Ep) q ■ Vdr 
4 qvp 


8uj 

q 
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where Ea{q,u}) is a complex function defined in terms of 
Za = {2 — aui')/q and Wa = ctA' jq as below. 


F„(g,A) = g[l-Z2] 


' [I + W^Zaf 

[1 - W2] [1 - Z2] 


(9) 


with vp = Sato/2h k, 10® to/s being the velocity of Dirac 
Fermions. The factor D{Ep) = gsgv\Ep\ /\/2TT{tivp)'^] is 
the DOS of the Dirac cones at E = Ep and the symmetry¬ 
breaking role of the electric current is manifested as the 
inner product of the wavevector q and the drift velocity 
of the electron (hole) gas, Vdr = sgn[Ep] vp {kshift/kp). 
The analytic expression for An(q,a;) given by Eq. (8) 
is the main result of this work and shown in Fig. 1. It 
conforms with the following principles: 
i) Real charge response. It satisfies the below condition 
which guarantees a non-imaginary charge response. 


An(-q,-a;) = [An(q,tc)]* (10) 


In principle, the occupation function of the drifting 
electron gas can be obtained via solving the Boltzmann 
transport equation (BTE)®^; however, in order to avoid 
the complexities of solving the BTE, we resort to the phe¬ 
nomenological shifted Fermi disk model which describes 
the nonequilibrium occupation function of a drifting elec¬ 
tron gas without the need for the details of the underlying 
scattering mechanisms®®40^ Pqj. ^ given shift of the Fermi 
disk from the Dirac point, kghift, the locus of the highest 
occupied states of the drifting electron gas with respect 
to the Dirac point, i.e. fc = 0, is given by, 

I ^1 - [^^^^sin0fc]2 - cos Ok] I (6) 

where kghift < kp is implied and the shift is assumed to 
be leftward, i.e. Eds || We thus limit ourselves to the 


ii) Causality. Since the integrand of AH has no poles 
in the upper half-plane of the complex frequency space, 
the real and imaginary parts of the analytic expression 
for An(q,w) are automatically correlated through the 
Kramers-Kronig (KK) relations, 

r [ = inAU{q,uj) (II) 

J UJ-UJ 

iii) The f-sum rule. The general f-sum rule for a bipartite 
tight-binding model® implies, 

J^AIl{q, uj)] OJ duj oe jE^k) Anp[E\k)](fk (12) 

where the right-hand side (RHS) represents the drift- 
induced modification to the kinetic energy of the tt elec¬ 
tron gas. The low-drift approximation for Anp{E,k) 
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FIG. 1. (Color online) Color-mapped values of (a) the real 
and (b) imaginary parts of An(q,a;) obtained from Eq. (8) 
and presented in units of {kahift/kF)D{EF). The positive and 
negative qx axes respectively correspond to the cases where 
q is anti-parallel {9q =0°) and parallel {Oq = 180°) to the 
drift velocity. The computed An(q,a;) values are corrected 
according to the Mermin’s approach^^ for a phenomenological 
scattering rate of h/r = 5meV (see Appendix B). 


given by Eq. (7) does not alter the kinetic energy of 
the electron gas, i.e. the RHS vanishes. Within the same 
level of approximation, the analytic expression for 3 [All] 
given by Eq. (8) is an even function of uj and yields a 
vanishing LHS, thereby satisfying the sum rule. 

It is worthy to note that the analytic expression for AH 
presented here contains only the terms that are linear in 
kshift/kp- The response for arbitrary drift velocity can 
be obtained numerically. In Fig. 2, the numerical result 
for the imaginary part of AH is shown for several (large) 
drift velocities at a fixed wave number q = 1.5kp in the 
direction of the drift. 

Let us finally note that within the framework of the 
shifted Fermi disk model, an exact analytic expression 
for the drift-induced intrasubband dynamical polariza¬ 
tion of two-dimensional electron gas (2DEG)^^ is obtain¬ 
able whose validity extends beyond the low-temperature 
and low-drift regime (see Appendix C), 

n*(q,a;) = n°(q,a; - ih/ml)[kshift ■ q]) (13) 

which suggests that the drifting 2DEG responds to the 



FIG. 2. (Golor online) The numeric-analytic comparison of 
A[An(q, ta)] for Dirac Fermions indicates that the analytic 
expression given by Eq. (8) becomes more accurate for smaller 
amounts of {kahijt/kF)- The numerical evaluation of AH is 
performed based on An_F values computed from Eq. (6), and 
the results are normalized by (kahift/kF) x D{Ef)- 


EM perturbation with a Doppler-shifted frequency. 


IV. DISCUSSION FOR DOPED SYSTEMS 
A. The static limit 


Within the low-drift and low-temperature regime, the 
analytic expression for the drift-induced modification to 
the intravalley static polarization of the tt electron gas in 
graphene can be obtained from Eq. (8) via setting a; = 0, 


An(q,a; 


0)^iD{Ep) 


q ■ Vdr 

Vp 



2 


q 

2k p 


]■ (14) 


One of the notable consequences of such a modification 
is the emergence of a drift-induced asymmetry in the 
Friedel oscillations (FO)^°. At far enough distances from 
the charged impurity atom, i.e. r>kp^, the modification 
to the FO in the presence of drift is described by, 

Auindir) ^ c sgn[Q] Vdr ■ f sm{2kpr) 

— 2 / 7 \2 

ris vf vf [KFry 

where ris = kp/n is the density of dopant electrons or 
holes, Q is the charge of the impurity atom, r is the 
in-plane position vector, c is the phase velocity of light 
in vacuum, af = e'^/A tteoHc (mks units) denotes the fine 
structure constant and kq is the background dielectric 
constant^°. Even though the drift-induced modihcation 
to the intervalley static polarization can be comparable 
to its intravalley counterpart, its contribution to the FO 
is negligible at far enough distances from the impurity 
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FIG. 3. (Color online) The computed (a) decay rate, (b) 
^piHpi ratio and (c) frequency of the TM-SPP modes of a 
suspended graphene channel, propagating parallel (red solid 
curves) and anti-parallel (blue dashed) to the drift velocity, 
are compared with their no-drift counterpart (black dotted 
curves). The agreement with the behavior predicted by 
the local approximation (green dash-dotted) in the q kp 
limit can be seen. The impact of the disorder-induced electron 
scattering is taken into account by the Mermin’s approach 
using a phenomenological scattering rate of fi/r = 5meV. 


atom, i.e. r>kp^. This is because the reiativeiy iarge 
vaiiey separation ieads to rapidiy-osciiiating terms in the 
summation yieiding the intervaiiey contribution. 


B. The local plasmonic limit 


Within the local limit, i.e. g ^ |£u| <C 1, the dynamical 
polarization of doped tt electron gas is given by^°, 

n{q,uj)^^D{Ep)(j-'^ (16) 


and the expression for AH given by Eq. (8) reduces to, 

(U) 

This translates into the drift-induced modihcation to the 
Drude weightgiven as follows, 

AD = (18) 

where Dq = e^\EF \ / (ttS^) is the bare Drude weight of the 
dopant electrons or holes. The local-limit expression for 
the dynamical polarization given by Eq. (16) yields a 
dependence for the TM-SPP frequency which is shown 
in Fig. 3(c). On the other hand, an additional acoustic 
branch with a linear dispersion, i.e. ui ~ Vsq^ emerges 
in a double-layer"*’'*^'*® or gated"*^ graphene system. The 
drift-induced change to the Drude weight can usually be 
neglected for the optical branch, i.e., w ~ y/g. For the 
acoustic branch, however, this correction might become 
observable if the sound-velocity, Vs , is comparable to the 
drift velocity Vdr- 


C. Modified plasmon dispersion and damping 


A two-dimensional (2D) electron gas confined between 
two dielectric media of dielectric constants ei and e 2 
supports transverse magnetic surface plasmon polariton 
(TM-SPP) modes whose dispersion is yielded by the so¬ 
lutions of the following equation*^’'*8,49^ 




2afhc 

q 


n(q,a;) 


(19) 


The retardation region is defined as the region in the 
(g,u;) plane near the dispersion of light^®’®®. Since the 
EM fields corresponding to the modes located in the 
retardation region are poorly-localized to the graphene 
sheet'*®d9^ 

we discuss the effects of drift out of this 
regime, i.e. q/uj ^ vp/c, where the LHS of Eq. (19) 
reduces to ei - 1 - 62 . To simplify our study and block the 
plasmon damping pathways such as the plasmon decay 
into the intrinsic {Huj ~ 195 meV) and extrinsic optical 
phonon modes within the frequency range of our inter¬ 
est (fiw < we assume the Fermi energy to be 

low enough (Eir^O.leE) and we limit our study to the 
case of a nonpolar substrate. Comparing the dispersion 
relation in the absence and presence of drift at a fixed g 
yields. 


n*(q,cu;,(g))=n“(g,w°;(g)) (20) 

which implies that the electric flux modifies the TM-SPP 
frequency. Within the low-drift regime, Eq. (20) yields 
the drift-induced modification to the TM-SPP frequency 
ujpi{q) and decay rate 7 p/(q), 


Aujpi{q) - iA7p/(q) 


An(g,a;°;(g)) 


OT(q,a;) 


doj 



( 21 ) 
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Within the local limit, Eq. (21) yields an expression for 
the drift-induced modification to the frequency of the 
TM-SPP modes that is valid within the local limit, 

Aujpi{q) ^ ■ Vdr ; g < fcj- (22) 

This Doppler-like frequency modification suggests that in 
the local limit the plasmonic charge density excitations 
are partially dragged by the drifting tt electron (hole) 
gas. As it is shown in Fig. 3c, the presence of drift 
causes the TM-SPP dispersion of tt electron gas to split 
into two branches each of which corresponds to the TM- 
SPP modes propagating parallel (P), i.e. 9q = tt, and 
anti-parallel (AP), i.e. 9q = 0, to the drift velocity. This 
splitting can be inferred from the following expression, 

Aujpi (q) - iAjpi {q) = [q ■ v^r] T(g) (23) 


-Unoccupied 

Occupied 


I I Unoccupied 
^ Occupied 



where 3?[T(g)],Q[T(g)] > 0. The computed TM-SPP 
branches that are presented in Fig. 3(c) suggest that 
the drift-induced frequency splitting, i.e. 2gndr5ft[Tr(g)], 
reaches its maximum near q = Qc with g^ denoting the 
onset of the Landau damping^°d6 ^ggg Appendix D). 

Moreover, according to Fig. 3(a), the TM-SPP mode 
propagating parallel (anti-parallel) to the drift velocity 
has a longer (shorter) lifetime comparing to the case in 
which the drift is absent. As it is implied by Fig. 3(b), 
such a drift-induced change in the propagation length, 
if measured in units of the mode wavelength, reaches its 
maximum for the modes with g Ri Qc. As a result, the 
short (few-wavelength) propagation length of plasmons, 
which is the main challenge of graphene plasmonics®^, 
can be increased via the application of a drain-source 
voltage. This conclusion, nevertheless, is based on the 
assumption that the device temperature is not affected 
by the presence of the drift. Otherwise, to assess the 
overall drift-induced change in the TM-SPP propagation 
length, the increase in the temperature-induced plasmon 
damping®^’®^’®"* resulting from the Joule heating of the 
current-carrying device should be taken into account. 
The higher decay rate for the AP branch, suggests the 
possibility of the application of DC current as a plas¬ 
monic brake to establish a one-way EM waveguide®®^®^. 

As was noted in Sec. Ill, only the terms linear in 
kshift/kp are retained in the analytic expression. Such 
an approximation produces unphysical results within a 
tiny neighborhood of the onset of Landau damping, i.e. 
Qc ± Sq and ujpi{qc) ± Sui where 6uj = vpSq oc Vdr ■ q- The 
slight dip in the decay rate curve and the exaggerated 
peak in the ojpi/^pi curve of the P plasmons, which are 
respectively presented by Fig. 3(a) and Fig. 3(b), are 
the inevitable consequences of such an approximation. 
To remedy this shortcoming, the terms proportional to 
{kshift/kp)^-'^ should be derived and included in the an¬ 
alytic expression. 


FIG. 4. (Color online) Depiction of the fe-space dynamics®® of 
the TT electron gas (the straight arrows along the cross-section 
of the Dirac cones with the ky = 0 plane) in an undoped 
graphene channel subjected to a drain-source voltage. The 
drifting occupants of the valence band whose group velocity 
has an opposite component to Eds migrate to the conduction 
band through the Dirac crossing. Ultimately, the migrants 
backscatter to the valence band after losing a quantum of 
energy ftfl via the emission of a phonon mode (wavy arrow). 


V. ZERO DOPING AND PLASMON GAIN 

Aside from the mode which was predicted when in¬ 
cluding the vertex corrections®®, undoped graphene does 
not support any TM-SPP modes at T = OK within the 
RPA.®"* Here, we numerically show that a high enough 
drain-source voltage along an undoped graphene channel 
enables the channel to support specific TM-SPP modes, 
even for the purely hypothetical case of T = OK. More 
importantly, the numerical results indicate the possibility 
of the emission of the low-energy {hw < 30 meV) and long 
wavelength plasmons. Similar proposals can be found in 
Refs. 33 and 34 and references therein. 

The crossing nature of the conduction and valence 
bands in graphene obligates the drifting electrons in the 
valence band to move up to the conduction band be¬ 
cause VfcA®(fc), which is semi-classically interpreted as 
the group velocity, is not well-defined at fe = 0 for a single 
Dirac cone. That is, a drifting electron passing through 
the neutrality point must travel to the other band. As 
is shown by Fig. 4, if the drain-source voltage is high 
enough, the migrant electrons lose a quantum of energy 
Ml by emitting a phonon mode and backscatter to the 
valence band. 

The intrinsic high-held transport properties of metallic 
single-wall carbon nanotubes (SWCNTs)®® is one piece 
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of evidence that proves that such a transport model is 
a physically relevant one. In each of the two bands, i.e. 
s = ±1, the nonequilibrium electronic occupation can be 
modeled by a 0fc-dependent Fermi energy as illustrated 
in Fig. 4 and described by the following expression, 

= jh^[\cos 9k\ - cos 6k] ;s = ±l. (24) 

Regarding the ID nature of the conical subbands of the 
metallic SWCNT,®° the model given by Eq. (24) may not 
be the best one to describe the similar situation within 
the 2D Dirac cones of graphene. However, any other 
possible model should also allow for nearly vertical, i.e. 
qa <C 1, electron-hole recombination processes, which is 
a guarantee for the generation of the long wavelength 
interband plasmons. 

Here, the dynamical polarization of the tt electron gas 
in a current-saturated and intrinsic graphene channel is 
approximated by feeding the response function with the 
nonequilibrium Fermi energy given by Eq. (24). In order 
to identify the well-localized TM-SPP modes in the (q.w) 
plane, we rely on the energy loss function S{q,uj), which 
is a measure of the spectral intensity of these modes®°. 


S{q,u}) = -3 



2afhc 

[ei + £ 2 ] q 


n(q,u;) 


(25) 


The computed energy loss, which is presented in Fig. 5, 
suggests that: i) the presence of electric current causes 
certain TM-SPP modes to emerge (which we refer to as 
the “drift-born” modes), ii) The electric current does 
not introduce any asymmetry to the response. The latter 
is due to the peculiar nonequilibrium occupation of the 
drifting tt electrons. The most important feature of the 
nonequilibrium response presented in Fig. 5 is the nega¬ 
tive energy loss of TM-SPP modes with qmin < 9 < qmax- 
Accordingly, we propose the use of the current-saturation 
in a nearly-undoped graphene channel as a mechanism for 
the amplification of THz plasmons. 

Since the electron-electron (e-e) interactions in the tt 
electron gas in graphene become significant for very low 
densities of dopant electrons,it is necessary to dis¬ 
cuss whether the validity of the results presented for 
the undoped case is challenged by the e-e interactions. 
Experimentally®^ as well as theoretically®^, it has been 
established that there is no gap opening due to chi¬ 
ral symmetry breaking even at the largest effective cou¬ 
pling constant present in suspended graphene (a^ = 2.2). 
However, strong Fermi velocity renormalisation takes 
place for low densities of the order Ug ^ 10^®cm“^ up 
to a factor of three.®®’®® To a first approximation, this ef¬ 
fect can be taken into account by using the renormalized 
Fermi velocity instead of the bare one (i.e., vp)- There¬ 
fore, the e-e interactions would not hinder the plasmon 
amplification mechanism proposed here but rather would 
modify the quantitative aspects of the numerical results 
we presented in the non-interacting picture (Fig. 5), and 
additional work is needed to clarify this issue. 
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FIG. 5. (Color online) Color-mapped values of the energy 
loss function, S{q,u}), of the tt electrons in a suspended and 
undoped graphene channel along which a high drain-source 
voltage is applied. The negative energy loss for the TM-SPP 
modes with qmin < 9 < qmax indicates the possibility of the 
amplification of these modes through the current saturation 
mechanism (AD « 149 meF). The negative and positive q^ 
axes respectively correspond to the P and AP cases. 


The electron and hole puddles, i.e. the spatial fluc¬ 
tuations in the Fermi energy, set a technical barrier 
to achieving a uniform neutrality along the graphene 
channel.®^ Fortunately, it has been shown experimen¬ 
tally that these charge puddles can be substantially re¬ 
duced on a hexagonal boron nitride (hBN) substrate,®®’®® 
thereby leaving some possibility for the plasmon am¬ 
plification mechanism proposed in this work. How¬ 
ever, theoretical investigations^®^^^ and experimental 
measurements^®’^'^ suggest that a proper crystallographic 
alignment of graphene with the hBN leads to the local 
breaking of the sublattice symmetry, thus opening a siz¬ 
able band gap at the Dirac point. There are two grounds 
on which it can be shown that the use of hBN substrate 
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does not necessarily induce any band gap: i) The sublat¬ 
tice symmetry in graphene can only be broken for specific 
relative rotation angles between the crystals, and this is 
why several works failed to detect such a gap.®®>®9d5,76 
Even if graphene is properly aligned with the hBN crystal 
so that the band gap emerges, placing an additional hBN 
crystal on top of graphene would kill the commensurate 
state and recover the sublattice symmetry7^d7 There¬ 
fore, such an unfavorable gap, which obstructs the pro¬ 
posed plasmon amplification mechanism, can be feasibly 
avoided by a proper encapsulation of graphene with hBN. 

Regarding the high temperature of a graphene channel 
within the current-saturation regime,^® it is necessary to 
incorporate the effects of temperature into the evaluation 
of the energy loss function, and as is shown in Fig. 5, the 
negative energy loss persists at high temperatures. 


VI. SUMMARY AND CONCLUSIONS 

We have discussed the dynamical response of a Dirac 
system subjected to a source-drain current. This was 
done by considering the nonequilibrium distribution 
function and feeding it into the well-known Lindhard 
function^®. By this, we were able to obtain closed-form 
expressions within the low-drift limit and analyzed the 
nonequilibrium response. Since the f-sum rule is obeyed 
in this limit, our approximation can be regarded as quasi¬ 
equilibrium. However, the sum rule does not hold any¬ 
more for the case of kshift ^ kp, i.e., for systems out of 
equilibrium. 

For doped systems, the asymmetric response of the 
drifting tt electron gas was discussed in the static and 
dynamical limit, especially commenting on the modified 
plasmon dispersion and damping rates. For a neutral sys¬ 
tem, where the external electric field does not induce an 
asymmetric response due to particle-hole symmetry, nu¬ 
merical results of the energy loss function were presented 
and a plasmon gain region found which persists even at 
high temperatures such as T = 600K. This result may 
be relevant and lead to potential applications based on 
ultra-clean encapsulated graphene samples. 

We finally discuss the limitations of the shifted Fermi 
sea model. Relaxation processes in graphene are known 
to be fast—in particular, e-e relaxation tends to equi¬ 
librate the system within femtoseconds. Our model is 
therefore mainly valid in the analytical limit kshift/kp 
1 which is reaffirmed by the fact that the sum-rule holds 
in this case of quasi-equilibrium. However, note that the 
Doppler-like transformation implied by Eq. (13) is in 
concordance with the experimental measurements on the 
drifting 2DEG system®^. This can be regarded as sup¬ 
porting evidence for our approximate treatment of the 
EM response of a driven electron gas. 

Several extensions are possible such as discussing 
the response of a gapped system, multi-layer systems 
and comparing our results with non-linear response 
functions.®^ 


ACKNOWLEDGMENTS 

We thank Daniel R. Mason and Guillermo Gomez San¬ 
tos for useful discussions. This work has been supported 
by the National Research Foundation of Korea (National 
Honor Scientist Program, 2010-0020414) and by the Min- 
isterio de Economia y Gompetitividad (FIS2013-48048-P, 
FIS2014-57432-P). 


Appendix A: Derivation of the drift-induced 
modification to the dynamical polarization 


For an n(a p)-doped graphene sample, the contribution 
of the occupied (empty) eigen-states of the conduction 
(valence) band to the dynamical polarization of the tt 
electron (hole) gas can be separated out as follows^°: 

n(q,w) - n^-=°(q,u;) = D{Ef) ^ t/„(q,u;) (Al) 

a—± 

At T = 0, the complex function Ua{q,uj) is given by 

(A2) 

, where Aa{k,q,oj) is defined as follows: 

A — i, f+,a{k,q) 

^^l3w+VF[k-a\k + q\]+il3Q+ ^ ^ 


However, in a case where the rotational symmetry around 
the Dirac point is broken (e.g., a dfc-dependent Fermi 
wavevector), the formalism given by Eqs. (Al), (A2) 
and (A3) does not satisfy the condition in Eq. (10). It 
is straightforward to derive the general formalism that 
holds for the non-symmetric case directly from Eq. (2). 
The appearance of jSq instead of q is the feature which 
distinguishes the general formalism from the old one. 


Aa = vpk 

0 =± 


__ 

fiui+VF [k—a \k+j3q\\+ij30'^ 


(A4) 


Hence, the drift-induced modification to the dynamical 
polarization is given by, 

An(q,c^) = rd0k /'^{a; + A'_}dk. (A5) 

Stt/cf Jo Jfep 

Within the low-drift regime, Eq. (6) reduces to 

fc^-fcF|l-^cosdfc| ;^«1- (A6) 

Regarding the small drift-induced perturbation to the 
Fermi wavevector, i.e. \kp — kpl kp, the fc-integral 
can be approximated as follows: 

rkF—kshift cos 6k 

I A^ dk = — kshift cos 0k (A7) 

J kp 







Substituting this result into Eq. (A5) yields 
D{Ep) kshift 


Jo 


where the integrand in Eq. (A8) is given by 
Ba{q,uj,ip) = -[A'^{k,q,u;)]^^f^^ cos {ip + Og) 


(AS) 


(A9) 


The real part of the integral in Eq. (A8) can be obtained 
via applying the following integral identity: 


^271- 


dip 


Iq I — p cos ip 


= 27r 


e[i-p^] 




(AlO) 


with p being a real number and 0 denotes the Heaviside’s 
step function. We then arrive at the following expression, 


5R[An(q,a;)] ^ A 


8w 

q 




a—dz 


(All) 





where the function G^{q,u]) reads 

|(:>(a)-2a)-g^| [((:}-2a;)^-g^] 

q^{Q^-q2)[(Q-2ar-<P] 


(A12) 


FIG. 6. (Color online) The analytic solution for An(q,cj) can 
be specified in each of the six regions in the {q, u) plane which 
are outlined by the straight lines of u) — q (solid), cD = 2 + g 
(dash-dotted), Co = 2 — q (dotted) and Co = q — 2 (dashed). 


and the coefficient is given by 

= -ba0[(w^ - q^){{ui - 2af - g^}] . (A13) 

The following integral identity leads us to the imaginary 
part of the integral in Eq. (A8): 


hence, to have a better presentation of the physical as¬ 
pects of this phenomenon, we rewrite the coefficient A in 
terms of the drift velocity 



N{ip)dp 
M{ip) + 



N{p) 

dM(ip) 

dip 


if^ipj 


(AM) 


with N{ip) and M{p) being two analytic functions within 
the range of [0, 27r], and pi ^2 are the duet roots of M{p). 
The resulting expression for 5[An] is as follows: 

S[An(q,w)] ^ ^^DlGi{q,Q) sgn[uo - a] (A15) 

® a=± 


where the real function Gjy,{q,Lj) is described by 

(A16) 

gV(9"-w2)[(w-2a)2-g2] 

and the coefficient Di reads as follows: 


Vdr = ±^fAnF[E^{k)]kd^k. (A19) 

J 

Plugging AnpiE, k) from Eq. (7), which corresponds to 
the case of kshift = —ex, into Eq. (A19) yields 

^ = sgn[Ep] . (A20) 

vp kp 

Rewriting the coefficient A in terms of the drift velocity 
from Eq. (A20) yields the prefactor appearing in Eq. 
(8). Fortunately, the expansive expressions for the real 
and imaginary parts of An(< 7 , uo) can be compacted into 
the complex function given by Eq. (8), which applies to 
complex frequency values. On the other hand, in each of 
the regions specified in Fig. 6, the real part of An(q, uo) 
can be expressed in terms of real functions: 


DI = -a 0[(g^ - (1)^) {(w - 2a)^ - q^]] (A17) 

Lastly, the coefficient A is given by 

A= \ sgn[Ep] D{Ep) cos Oq . (A18) 

4 kp 

The sgn[Ep] factor indicates that the drift-induced 
asymmetricity in the Ep >0 (E^’ < 0) case is solely dic¬ 
tated by the drift velocity, Vdr, of the electron (hole) gas; 


K[An] 


8Auj A^ 

——I—. < 

q x/|w2-g2| 


0 

lA 

H+ + H^ 

IB 

-H+ 

2A 

+i7- 

2B 

H_-H+ 

3A 

H_-H+ 

W 


(A21) 


Likewise, the imaginary part of AH)*/, uo) is given by 
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An^ 

3[An]^—^J=X< 

^/W^\ 


Ha{q,Q) = { 1 + 


auj 


2 —aw 



lA 

0 


IB 

H_ 

2A 

H+ 


2B 

0 

3A 

0 

3B 

given by 


\ 

1- 

2 — auj 

i\ 

_ q 


(A22) 


(A23) 


Having evaluated AH in the {q, w > 0) quarter-plane, we 
can evaluate AH in the (g, —w < 0) quarter-plane by ex¬ 
ploiting the fact that the real and imaginary parts of AH 
are odd and even functions of w, respectively. Note that 
this is just the opposite symmetry relation as compared 
to the equilibrium case. 


Appendix B: The Mermin’s approach 


Regarding the complexity of the electron scattering 
mechanisms, their impact on the response function can 
be approximately taken into account by replacing w in 
the response function of the collisionless electron gas by 
w -f iT~^ with T~^ being the phenomenological electron 
scattering rate which is related to the mobility of the 
graphene sample ^ via the following relation:®^ 


evl 


(Bl) 


Such an imprecise scheme, however, fails to conserve the 
local electron number. The following correction formula 
removes such a defect for the case of the intrasubband 
longitudinal response function of the 2DEG^^ 


n^(q,w) = 


n(q,w-|-iT 


1 - 


1 — iujT 


1 - 


n(q,w-|-iT 


n(q,0) 


(B2) 


where n(q,w) and ni-(q,w) denote the collisionless and 
the corrected dynamical polarization, respectively. Such 
a correction scheme has been shown to be applicable 
to the intraband dynamical polarization of graphene®®; 
however, to our knowledge, there is no literature in which 
the application of Mermin’s approach to the interband 
dynamical polarization of graphene is rigorously justified. 
Nonetheless, we relegate the clarification of this matter to 
the future works and follow the general trend of applying 
the Eq. (B2) to the case of Dirac Fermions "*’®’®^“®®. 

The dynamical polarization of the tt electron gas in 
a doped graphene sample at T = QK is given by the 
following complex function:®^ 




8v^ 





(B3) 


where = (2 — aCj')/q and the complex function G°^{z) 
is defined as below: 

G°‘{z) = z\J\ — z'^ + ai\\i[z + \/ z — 1\/ z + 1] (B4) 

The effects of the disorder-induced electron scattering on 
the dynamical polarization of non-drifting tt electron gas 
in graphene can be taken into account by feeding the 
Il{q,uj + iT~^) values from Eq. (B3) into Eq. (B2). To 
obtain H* (q, w), we have computed An(q, w-f-fr”®) using 
Eq. (8), added it to \i{q,uj + iT~^) given by Eq. (B3), 
and fed their sum into Eq. (B2). The H* — n° val¬ 
ues computed for a phenomenological scattering rate of 
h/r = bmeV are presented in Fig. 1. For a Fermi en¬ 
ergy of 100 meV, this r-value corresponds to a sample 
mobility oi 10^ as suggested by Eq. (Bl). 


Appendix C: The case of the 2DEG 


Feeding the integral in Eq. (2) with a single parabolic 
band, i.e. E{k) = h^k^/2ml along with excluding 
and fs,s' (fe, q) yields the intrasubband dynamical polar¬ 
ization of the two-dimensional electron gas (2DEG)'*^. 
Thanks to the absence of interband transitions and the 
parabolic energy dispersion, a change of the integration 
variable according to k' = k~kghift yields an explicit re¬ 
lation between n*(q,w) and n°(q,w) given by Eq. (13). 
In spite of the fact that the transformation given by Eq. 
(13) provides the best possible approximation within the 
framework of the shifted Fermi disk model, we present 
the 2DEG counterpart of the expression given by Eq. 
(8) in order to have a comparative picture: 


An(q,w) 


q ■ Vdr 
qVF 


D{Ef) 

q 



a 


\ r 

_Lj' — aq^ _ 


(Cl) 


with D{EF) = gsml/2Trh'^ being the DOS of the parabolic 
band and Vdr = VF[kshift/kF] is the drift velocity. 


Appendix D: The onset of the Landau damping in 
the absence of drift 

In the absence of drift, the TM-SPP dispersion can be 
implicitly obtained via plugging the analytic expression 
for the dynamical polarization given by Eq. (B3) into 
Eq. (19). The onset of the Landau damping qc is the 
g-value where the TM-SPP dispersion curve and the line 
uj = 2 — q meet which is yielded by the below equation: 

1^1= G> (1 -1) = 1 + AAiliA (Di) 

16 Vl - gc \qc ) 2gsgvaf{c/vF) 

where the function G>(a;) is given by 

G> (x) = x\/— 1 — In [a; -f \/x"^ — 1] . 


W[q,oj) = D{EF) 


(D2) 
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